package de.latzko.pattern.monteCarlo

import java.util.Random
import java.util.Random._
import de.latzko.pattern.math._
import de.latzko.pattern.math.NearestNeighbour
import scalala.tensor.dense.DenseVector._
import scalala.tensor.dense.DenseVector
import scalala.Scalala._
import scalala.Scalala

/**
 * 
 * @author tom
 * vectors stores the data
 *	q is the number of spins
 *	k nearest neighbors 
 */
class SwendsenWang(val vectors: Array[DenseVector], val q: Int, val k: Int) extends BaseMonteCarlo
  with NearestNeighbour with Mesure  {

	// variables 
  var bonds = new Array[Byte]((vectors.size * vectors.size - 1) / 2);
  var label = new Array[Int](vectors.size);
  var spins = new Array[Int](vectors.size);
  val n = vectors.size;
  var rand = new Random();
  var vecs = vectors.map(x => (knnPlusDist(x, vectors, k)))
  var meanAvr: Double = math.pow((minAverageDistance + maxAverageDistance) / 2, 2)

  val getBondsIj = (i: Int, j: Int) => bonds(i + (j * (j - 1)) / 2)
  val calcJij = (dist: Double, avr: Double, k: Int) => { 1.0 / k * math.exp(-dist / avr) };
  vecs = vecs.map(x => x.map(y => { (calcJij(y._1 * y._1, meanAvr, k), y._2) }))

  /**
   * Start the methode
   * @param mcSteps compute mcSteps 
   * @param T 		at temperature T
   * @return
   */
  def runMethode(mcSteps: Int, T: Double): Double = {
    init()
    val thermSteps = mcSteps
    for (i <- 1 to thermSteps)
      monteCarloStep(T);
    for (i <- 1 to mcSteps) {
      monteCarloStep(T)
      mesure(spins, n, q);
    }
    val moments = computeAverages()

    return (n / T) * (moments._2 - (moments._1 * moments._1))

  }

  /**
   * compute one montecarlo step
   * @param T temperature
   */
  def monteCarloStep(T: Double) = {
    //  var t = System.currentTimeMillis();
    meltOrFrozen(T);

    //   println("melt " + (System.currentTimeMillis() - t))
    //  t = System.currentTimeMillis()
    labelSwCluster();

    //  println("label " + (System.currentTimeMillis() - t))
    // t = System.currentTimeMillis()
    flipSpins();

    //  println("flip " + (System.currentTimeMillis() - t))
   // scatterPlot(vectors, label)
    label = label.map(x => 0)
  }

  /**
   * bring system in start 
   */
  override def init(): Unit = {
    super.init()
    spins = spins.map(x => rand.nextInt(q))

  }

  /**
   * Create a field with conntected nodes
   * @param T is the temperature 
   */
  private def meltOrFrozen(T: Double): Unit = {
    bonds = bonds.map(x => 0.asInstanceOf[Byte]);

    for (i <- 0 until n)
      vecs(i).foreach(x => freezeBonds(i, x._2, x._1, T))

  }

  /**
   * flip the spins of a hole cluster via random spin
   */
  private def flipSpins(): Unit = {
    val maxLabel = label.reduceRight(_ max _)
    var randArr = new Array[Int](maxLabel)
    randArr = randArr.map(x => rand.nextInt(q))
    var i = -1;
    label.foreach(x => { i += 1; if (x > 0) spins(i) = randArr(x - 1) })
  }

  /**
   * decide that i and j are frozen together
   * @param i
   * @param j
   * @param Jij
   * @param T
   */
  def freezeBonds(i: Int, j: Int, Jij: Double, T: Double): Unit = {

    val pos: Int = i + j * (j - 1) / 2;

    if (bonds(pos) == 0) {
      var res: Byte = 1;
      if (spins(i) == spins(j)) {
        val freezeProb = 1 - math.exp(-Jij / T);
        if (rand.nextDouble() < freezeProb)
          res = 2;
      }

      if (i > j) {
        bonds(pos) = res;
      } else {
        bonds(pos) = res;
      }
    }

  }

  /**
   * find the cluster that are frozen together
   */
  def labelSwCluster(): Unit = {

    var newLabel = 1
    for (i <- 0 until n - 1) {
      if (label(i) == 0) {
        val connected = growSwLabel(i, newLabel);
        if (connected)
          newLabel += 1;
        else
          label(i) = 0

      }
    }
  }

  /**
   * grow the cluster 
   * @param i
   * @param newLabel
   * @return
   */
  def growSwLabel(i: Int, newLabel: Int): Boolean = {
    var connected = false;
    if (label(i) == 0) {
      label(i) = newLabel
      for (j <- i + 1 until n) {

        if (getBondsIj(i, j) == 2) {
          if (label(j) == 0)
            growSwLabel(j, newLabel);
          connected = true
        }
      }
    }
    return connected
  }

}